Correlation Energy Divergences in Metallic Systems 
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We numerically examine divergences of the total energy in metallic systems of approximate many- 
body theories using Hartree-Fock as a reference, including perturbative (Mpller-Plesset, MP), cou- 
pled cluster (CC) and configuration interaction (CI) approaches. Controlling for finite size effects 
and basis set incompleteness error by comparison with energies from the random phase approxima- 
tion (RPA), we suggest convincingly that non-perturbative coupled cluster theories are acceptable 
for modelling electronic interactions in metals whilst perturbative coupled cluster theories are not. 
Data are provided from the RPA with which it is possible to test other approximate finite-basis 
methods for divergences with only modest computational cost. 
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Introduction. - Although it has long been known that 
finite order perturbation theory predicts divergent terms 
in the energy for metals, as analysed by Gell-Mann and 
Brucckner in 19570, it has taken fifty years for this 
result to be reproduced numerically for a real metallic 
system Q. This is in part because simulating the many- 
body wavefunction of solid state materials starting from 
Hartrcc-Fock theory is incredibly computationally ex- 
pensive, often scaling polynomially with a variety of sim- 
ulation parameters including the size of the simulation 
cell. In spite of the difficulties in applying post- Hartree- 
Fock methods to solids, there is an increasing body of au- 
thors attempting to do so with a variety of wavefunction 
ansatzes on the so-called 'hierarchies of quantum chem- 
istry' Q 0, 0-[IH ■ Their aim is to achieve highly accurate 
correlation energies for the field of materials science by 
leveraging the systematic improvability of these approx- 
imate solutions to the Schrodingcr equation 1Q, [17 1 . 

An open question of growing importance surrounding 
this field is to directly address which methods are ap- 
propriate and which are not for the study of metallic 
systems [lij]. Approximations and divergences need to be 
understood so that needless effort is not expended inves- 
tigating methods which will ultimately fail. Although 
it would in principle be possible to pursue this ques- 
tion with analytical theory, the plurality of diagrams 
and the lack of closed solutions makes this attempt in- 
tractable, especially for higher orders of coupled cluster 
theory. So far most analytical results have been achieved 
by neglecting the exchange term in the Hartree-Fock ref- 
erence which, in perturbative theories, physically corre- 
sponds to starting with a non-interacting reference. This 



greatly simplifies the expression for the denominator in 
M0ller-Plesset perturbation th eoryll9l and allows analyt- 
ical divergences to be found [lL l20l 12 1| . Strictly speaking, 
however, these results consider only a lower bound to 
the true correlation energy of the respective method us- 
ing the Hartree-Fock reference. On the other hand, one 
might think that a closing band gap is sufficient to guar- 
antee divergences in a perturbative approach, however, 
it is well-known that some terms converge |22f . Further- 
more, especially in real systems, numerical divergences 
can be hard to find and distinguish from convergences. 
To this end we aim to provide here a simple, novel and 
robust methodology to test for the numerical convergence 
of approximate methods in metals using the finite basis 
set simulation-cell electron gas [23H25l |. 

We focus in this study on testing a number of famous 
and commonly-used methods at the heart of quantum 
chemistry. These are coupled cluster theory with full am- 
plitude equations [1^] (l7ll2^ - l30| . M0ller-Plesset perturba- 
tion theory [Tijl and truncated configuration interaction 
(singles and doubles) (31. 32 1 . These methods have been 
extremely successful in treating correlation in molecular 
systems and they are being increasingly applied to solid- 
state systems; extension to metallic systems is an ongoing 
aim in this community. We also draw comparison with 
the random phase approximation [33j that has already 
received a gr eat deal of attention from the solid-state 
physicists 3414361 ] . These are only presented as examples, 



the methodology here is extendible to any wavefunction 
method which can be formulated in a finite canonical 
Hartree-Fock basis 37-391. 
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Model. - The homogenous electron gas (uniform elec- 
tron gas) consists of N electrons in a box of length L 
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with a two-electron Ewald interaction v a a 



const. 



(1) 



In the thermodynamic limit, found as the particle num- 
ber tends to infinity (TV — > 00) with the density held 
constant, it is possible to solve the above Hamiltonian in 
the Hartrcc-Fock Approximation with plane waves. This 
yields an analytic expression for the dispersion relation, 
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and x = k/kp, where this term is due to exchange |42j. 
It is this term which produces a band structure with a 
zero in the density of states at the Fermi energy, and also 
complicates analytical derivations. 

In setting out to find the behaviour of approximate 
theories to obtain the correlation energy (i.e. the to- 
tal energy with Hartrcc-Fock energy as a starting-point), 
it is typical to start with a finite simulation-cell model 
of TV electrons, and carefull y appr oach the thermody- 
namic limit by extrapolation|4ll |43|. However, in quan- 
tum chemical techniques, we must also make do with a 
finite one-particle basis set. Very little has been studied 
about the relatively simple properties of the HEG repre- 
sented in a finite basis, at least in part since the study 
of the electron gas has been pushed forward so much by 
continuum quantum Monte Carlo techniques that work 
at the complete basis set limit (44- 46 1. The difficulty of 
investigating the properties of these approximate theories 
in the thermodynamic limit is hampered by this require- 
ment of a finite basis set, in this case of M plane waves 
spinorbitals defined by a kinetic energy cutoff \k\. In 
principle, the complete basis set limit k c — > 00 and ther- 
modynamic limit TV — > 00 must be found, which is pro- 
hibitively costly given the scaling of even approximate 
quantum chemical theories. 

The most obvious way to make progress towards these 
limits is to take the k c — > 00 limit, to solve the TV-electron 
Hamiltonian at the complete basis set limit, and then the 
TV — >• 00 limit can be found latterly. However, in this 
study we p ropose to take the TV — > 00 limit first for a fi- 
nite fc c [47|. We first outline how to show the well-known 
divergence in the MP2 energy using finite-M, finite-TV 
calculations and then generalise this approach to demon- 
strate limiting behaviours in other theories. 

In the thermodynamic limit, the 3D electron gas mo- 
mentum space is a Fermi sphere containing all momenta 
less than the Fermi wavevector k < kp. The virtual man- 
ifold of free-electron plane waves surround this sphere 
stretching out to infinite k and infinitesimal spacing. In 



simulating this with finite systems, the finite electron 
number controls the spacing of the allowed plane-wave 
wavevectors and a finite basis set is produced by a cut- 
off, \k 2 < \k 2 c , which controls the extent to which the 
plane waves stretch out into fc-space around the T-point. 
Allowing TV — > 00 for a 'finite-basis' electron gas is well- 
defined for a given 7 = k c /kp and amounts to represent- 
ing the space inside k < k c with an ever finer grid. 

We can then simulate a series of TV-electron gases with 
M basis functions where M = 7 3 TV to within finite-size 
effects jijj]. For a given 7, as the TV — > 00 limit is taken, 
the band gap closes because the grid spacing in the re- 
gion around the Fermi surface becomes smaller, and the 
zero-momentum (q = 0) excitations that cause the di- 
vergences in for example MP 2 theory are increasingly 
well-represented in a size-consistent fashion [49j. 

Results. - MP2 correlation energies are presented in 
Fig. [T] for sets of finite- TV, finite-M electron gases con- 
structed in this fashion. This conclusively demonstrates 
that this approach recovers the expected divergence and 
physical behaviour from this method which is insensitive 
to our choice of 7. In contrast, the differences between 
the different values of 7 converge to give a finite energy 
in the thermodynamic limit (TDL) represented by the 
lines for different 7 being parallel to one another within 
extrapolation error. Technically, this is just the basis set 
incompleteness error, which should be recovered as I/7 3 . 
To further validate this as a physical approach accurately 
capturing the TDL, we compare this with the finite-basis 
electron gas energies from identically constructed RPA 
calculations, which show a convergent behaviour (with a 
finite-size error as ~ TV" 1 51) as anticipated. All RPA 
results in this work are calculated using a HF reference. 

We further note that the band gaps of the simulation- 
cell electron gas in these simulations, at r s = 1.0 (rela- 
tively high density), arc still fairly large (54.8 — 6.6 eV) 
due to remaining finite-size errors even though we are 
simulating electron numbers between TV = 14 — 1030. 
To examine the insensitivity of the divergent behaviour 
with respect to a wider range of band gaps, we have 
compared simulations of several densities (by changing 
the r s parameter) over the same range of electron num- 
bers in Fig. [2j As the density is lowered and r s is 
raised, the range of typical metallic densities of metals 
is traversed [13, and although different contributions of 
kinetic and exchange energy are present in the Hartrcc- 
Fock energy, the divergence is remarkably unchanged. 
Furthermore, we now see that the divergence persists at 
a very large range of simulation-cell band gaps. This un- 
derlines a main conclusion of this study, that in spite of 
simulating what appear to be insulating finite systems, 
the divergence in the MP2 energy is still visible and the 
energy is without even qualitative physical relevance. 

Having now numerically demonstrated the well- 
accepted behaviour for the MP2 energy, we turn our at- 
tention towards other approximate methods for which 
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FIG. 1: (Colour online) MP2 energies for a variety of 
finite-basis simulation-cell electron gases with electron 

numbers N = 14 — 1030 corresponding to closed-shell 
configurations of a simple cubic reciprocal-space lattice 
and density r s = 1.0 a.u. These all diverge as 1/E gap as 
the thermodynamic limit (TDL) is approached and the 

three different basis set sizes are parallel to within 
errors from the fit. For a single basis set size, the RPA 
correlation energy is also shown, which converges in the 

TDL. 



the behaviour on approach to the TDL is unverified - 
coupled cluster singles and doubles theory (CCSD) and 
the addition of perturbative triples (CCSD(T)). There 
has been surprisingly little literature concerning this hi- 
erarchy of methods as applied to solids, in spite of the 
wealth of applications they have received in the molec- 
ular quantum chemistry community. Even though there 
has been some discussion of CCSD with approximate am- 
plitude equations, these more closely resemble the RPA 
equations [2l[ I51I453J. As such, to the best of our knowl- 
edge, the question of whether CCSD and CCSD(T) di- 
verge in the TDL for metallic systems has not yet been 
conclusively addressed and requires further investigation. 
We note in passing that CCSD and coupled cluster dou- 
bles theory (CCD) are equivalent for the homogeneous 
electron gas due to the complete absence of symmetry- 
allowed single-excitations in its many-body expansion. 

Due to the relatively expensive scaling of such meth- 
ods, simulations of a N = 1030 electron gas with current 
fullv-periodic codes [54^ are prohibitively expensive. How- 
ever, we have found that further reduction in finite-size 
effects can be achieved by taking the difference between 
the CCSD or CCSD(T) energy with the RPA energy and 
in this difference the limiting behaviour is more clear. 
We have also taken advantage of other simulation-cell 
lattices (face-centered cubic and body-centered cubic) to 
provide more closed-shell configurations. Taking energy 
differences in this way allows us to clearly demonstrate, 
in Fig. [3J that the CCSD energy is strongly convergent, 



FIG. 2: (Colour online) MP2 energies for a single basis 
set size 7 = \[2 for N = 168 — 1030 at a variety of 
densities showing that the divergent behaviour is 
insensitive to the density parameter. To compare on the 
same scale, the energies in both axes have been 
multiplied through by r s , such that for example the 
true band gaps of the r s = 10.0 a.u. simulations span 
the energy range ~ 0.5 — 1.0 eV. 



whereas the CCSD(T) diverges at the same speed as the 
MP2 energy [H3. 

For a better understanding of the above results we may 
consider the dominant Goldstone diagrams included in 
the respective methods. Fig. [4] shows the diagrams re- 
sponsible for the divergences in MP2 and CCSD(T) as 
well as an diagram contributing to the CCSD energy. 
The horizontal (double) wavy lines correspond to the 
(screened) Coulomb interactions. The vertical lines de- 
note occupied and unoccupied orbitals. Algebraic ex- 
pressions given beneath these diagrams correspond to 
the HEG assuming a point-like Fermi sphere[56j. Evi- 
dently, the unscreened Coulomb interactions are singu- 
lar for zero momentum-transfer vectors. In combination 
with the closing Hartree-Fock (HF) gap, this is the rea- 
son for the divergence in the MP2 energy. In contrast, for 
CCSD, our numerical results suggest that the effectively 
screened Coulomb interaction is non-singular and leads to 
a finite CCSD correlation energy in metals. The effective 
non-singular Coulomb interaction is obtained by a sum- 
mation of infinitely many diagrams, which in the RPA 
correspond to the so-called ring diagrams and in CCSD 
is carried out by solving the amplitude equations (57j|. 
CCSD(T), however, includes a fourth-order diagram that 
corresponds algebraically to an expression that involves 
an integration over two unscreened Coulomb interactions 
and a closing HF gap. This expression closely resembles 
the MP2 expression. As such, the CCSD(T) correlation 
energy diverges at the same rate as the MP2, albeit with 
a smaller prcfactor as is corroborated by our numerical 
findings. 
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FIG. 3: (Colour online) Energy differences between the RPA correlation energy (r s = 1.0 a.u., 7 = v2) show that 
MP2 and CCSD(T) are divergent as the band gap closes (on approach to the thermodynamic limit). CCSD is 
convergent to a constant energy offset with respect to RPA which is only serendipitously close to zero for this r s and 
7. CISD gradually yields zero correlation energy per electron, due to a lack of size-extensivity. 
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FIG. 4: Goldstone diagrams present in MP2, CCSD and CCSD(T). Below the diagrams the respective algebraic 
expressions are given and approximated for the uniform electron gas assuming a Fermi point. k and a, &, c, d 
refer to occupied and unoccupied orbitals, respectively. Aeij a b and e g refer to the differences between the occupied 
and unoccupied one-electron energies. V c f{ denotes a non-singular screened Coulomb interaction. The exact form of 
V e ff depends on the form of the coupled cluster amplitude equations. In the RPA VJ.fr is given by e^ 1 (q)/q 2 , where 

e(q) is the dielectric matrix. 



For completeness, we also demonstrate the effect of this 
difference when the test method is not size-extensive. As 
an example of this, we used configuration interaction sin- 
gles and doubles (CISD, and equivalent to CID for the 
HEG), obtained by usin g a truncation of a CI quantum 
Monte Carlo calculation 58 - 6lJ , and this shows a ten- 
dency to retrieve an ever-lower fraction of the correlation 
energy as the electron number is raised (Fig. [3]). A simple 
rationalisation of this is as follows. When two A^-particlc 
systems are individually treated with CID, they achieve 



a better treatment of correlation than the combined 2N- 
particle system since in the individual systems, the con- 
sideration of quadruple excitations in the larger system 
is made absent 



31 



Due to this, correlation energy is 
not recovered consistently as N grows but, per particle, 
falls as some inverse polynomial. In particular, correla- 
tion energy is lost as N~z in the case of non- interacting 
electron pairs [l6, 62 1. 

Concluding remarks. - In summary, we have shown 
that a judicious choice of finite-size and finite-electron 
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number homogeneous electron gas models can be used to 
demonstrate the limiting behaviour of the correlation en- 
ergy in approximate many-body theories for metallic sys- 
tems with modest computational cost. By comparing to 
RPA correlation energies we control for basis set incom- 
pleteness and finite-size errors. As a first application of 
the outlined methodology, we have verified the divergence 
of MP2 energies in metals. Moreover, we have shown 
that CCSD(T) also exhibits a divergent behaviour in the 
HEG, by virtue of the perturbative (T) correction to the 
CCSD correlation energy. In contrast, CCSD, due to 
an "effectively screened" Coulomb interaction, predicts 
converging correlation energies in metals. Our findings 
strongly suggest that the cutting-edge accuracy of high- 
level perturbative quantum chemistry methods such as 
CCSD(T) will only be transferable from molecular quan- 
tum chemistry to the study of metallic solids, if novel 
theories are introduced that lift the divergent behaviour. 
This work provides a simple but stringent test for such 
novel many-body approximations that are too complex 
for analytical theory, and we have provided data such 
that other authors may conduct the same analysis [63 1 . 
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^1 7P-P"V"1"PT1 QI'tT'P 

energies 


Hartree-Fock Theory (HF) [HH 


N 3 


1st order 


m 


Yes 


O J J TV /T 11 T>1 j_ rr-ii /" Ti TT"> o \ |-i nil 

second-order Moller-Plesset lheory (MP 2) |19l| 


N 5 


2nd order 


Ul lol [7nl 


Yes 


Random phase approximation (RPA) [33|] 


N 4 


1st order 


MM 


Yes 


Coupled cluster doubles with approximated 
amplitude equations 


- 


approx. 3rd 
order 


MM 


Yes 


Coupled cluster singles and doubles theory (CCSD) 

[HQ 


N 6 


3rd order 




Yes 


Coupled cluster singles and doubles theory with 
perturbative triples (CCSD(T))[27[ 


N 7 


4th order 


BOJ, and[ll 


Yes 


Configuration interaction truncated at single and 
double excitations (CISD)[H 


N 6 






No 



TABLE I: (Supplemental Material.) Pedagogical summary of the acronyms, references and scaling of the canonical 
formulation of the employed many-body theories. We note that the scaling of these methods might be naturally 
lower for the HEG due to momentum symmetry constraints on the allowed excitations, but this is not considered 
here. References to solid-state applications are intended as examples rather than an exhaustive list. 




FIG. 5: (Colour online, supplemental material.) Scale diagram illustrating various finite-basis TDL electron gases at 
a variety of 7. The dark shaded circle indicates the location of the fermi sphere, and hence wavevectors of occupied 
orbitals, and the light section indicates wavevectors of virtual orbitals. 
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JjjlcCLIUIl IlUIIUJfcJI JV 


JTlcLIlc WdVc hJJillUI UlLcLlb Ivl 


fir UdllQ gcLJJ -Cvgap ^e V J 


V_>UIIclcLLlUIl eiieigy -C^corr. V * / 

electron) 


14 


38 


54.752 


-0.541 


38 


66 


32.521 


-0.260 


54 


162 


28.776 


-0.593 


66 


186 


24.357 


-0.675 


114 


294 


18.911 


-0.770 


162 


358 


24.614 


-0.565 


186 


514 


14.773 


-0.787 


246 


682 


12.642 


-0.830 


294 


778 


11.932 


-0.846 


342 


922 


11.787 


-0.827 


358 


970 


11.428 


-0.836 


406 


1174 


10.297 


-0.919 


514 


1502 


8.6289 


-0.915 


610 


1694 


8.5071 


-0.922 


682 


1850 


8.4075 


-0.904 


730 


2042 


8.2846 


-0.924 


778 


2090 


7.6773 


-0.903 


874 


2378 


7.9345 


-0.888 


922 


2474 


10.002 


-0.874 



TABLE II: (Supplemental Material.) RPA correlation energies and band gaps for the systems electron gas models 
studied here (r s = 1.0 a.u., 7 = \/2). Other lattice types are available on request by email. 



MP2 CCD 




MP2 CCD 




